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ABSTRACT 

By numerically evaluating five- dimensional phase-space integrals, we have calculated the 
small effect of finite nucleon mass on the weak-interaction rates that interconvert protons 
and neutrons in the early Universe. We have modified the standard code for primordial 
nucleosynthesis to include these corrections and find a small, systematic increase in the 4 He 
yield, AY/Y ~ (0.47 — 0.50)%, depending slightly on the value of the baryon-to-photon ratio 
rj. The fractional changes in the abundances of D, 3 He, and 7 Li range from 0.06% to 2.8% 
for 10- 11 <r]< 10~ 8 . 



1 Introduction 



Primordial nucleosynthesis is one of the cornerstones of the hot big-bang cosmology. The 
consistency between its predictions for the abundances of D, 3 He, 4 He and 7 Li and their 
inferred primordial abundances provides its earliest test. Further, big-bang nucleosynthesis 
has been used to obtain the best determination of the baryon density |], |[ || and to test 
particle-physics theories, e.g., the stringent limit to the number of light neutrino species 

Scrutiny of primordial nucleosynthesis, both on the theoretical side and on the obser- 
vational side, has been constant: Reaction rates have been updated and their uncertain- 
ties quantified [/], ||, 0; finite-temperature corrections have been taken into account [[10 



the effect of inhomogeneities in the baryon density explored |TTJ; the slight effect of the 
heating of neutrinos by e ± annihilations has been computed O, O]; the primordial abun- 



dance of 7 Li has been put on a firm basis and its destruction in stars has been studied 
T4], ITR 111, [17], [18], [T91 EH, |2H|; the production and destruction of 3 He and the destruction of 



D have been studied 0, ^2], |23], |2^] ; and astrophysicists now argue about the third significant 
figure in the primordial 4 He abundance [03, ^Bl 

A measure of the progress in this endeavour is provided by the shrinking of the "concor- 
dance region" of parameter space. The predicted and measured primordial abundances agree 
provided: the baryon-to-photon ratio lies in the narrow interval 2 x 1(T 10 <i]<7x 1(T 10 
and the equivalent number of light neutrino species N v < 4 j|, ||. The shrinking of the con- 
cordance interval motivates the study of smaller and smaller effects. In particular, once the 
primeval deuterium abundance is accurately determined in high-redshift hydrogen clouds, 
the baryon-to-photon ratio will be pegged to an accuracy of order 10%, and recent progress 
suggests that this will happen sooner rather than later p8| . In turn, this will reduce the 
uncertainty in the predicted 4 He abundance due to the baryon density to AY ~ ±0.005. 

The weak-interaction rates govern the neutron-to-proton ratio and thereby are crucial to 
the outcome of nucleosynthesis; e.g., the mass fraction of 4 He produced is roughly twice the 
neutron fraction at the time nucleosynthesis commences (T ~ 0.08 MeV). In the standard 
code [p9| , |30|| these rates are computed in the infinite-nucleon-mass limit because this sim- 
plifies the expressions for the rates to one- dimensional integrals. The finite-mass corrections 
involve terms of order m e /m^r, T/m^, and Q/m^, which are all of the order of 0.1%. Here 
m e is the electron mass, is the nucleon mass, T ~ O(MeV) is the temperature during 
the epoch of nucleosynthesis, and Q = m n — m p = 1.293 MeV is the neutron-proton mass 
difference. As it turns out, the corrections to the rates are actually of order a few percent, 
so that the change in 4 He abundance is expected to be a few parts in the third significant 
figure (AY/Y ~ —0.55 rate/rate). 

Because the third significant figure of the primordial 4 He abundance is now relevant, 
we set out to calculate the finite-nucleon-mass corrections to the weak-interaction rates 
by numerically integrating the exact expressions for the rates. This involved accurately 
(< 0.1 %) evaluating five-dimensional rate integrals. We incorporated our results into the 
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Figure 1: Weak rates as a function of temperature (infinite-nucleon-mass limit). Note, 
freeze-out of the n/p ratio occurs at T ~ 0.8 MeV and nucleosynthesis begins in earnest at 
T ~ 0.1 MeV. 

standard nucleosynthesis code and found that AY/Y « (0.47 — 0.50)%, depending on the 
value of rj. In the next Section we discuss the weak rates, and in the following Section 
finite-mass corrections to the weak rates. In Sec. IV we discuss our changes to the standard 
nucleosynthesis code, and we finish with a discussion of our results for the changes in the 
yields of 4 He and the other light elements. 



2 Weak- interact ion Rates 

The weak interactions that interconvert neutrons and protons, n p + e + u, n + e <->■ p + v, 
and n + v <-> p + e, play a crucial role as they govern the neutron fraction, and the neutron 
fraction ultimately determines the amount of nucleosynthesis that takes place. (Here and 
throughout we use e to indicate electron or positron, and v to indicate electron neutrino 
or antineutrino; the appropriate particle or antiparticle designation follows from charge and 
lepton- number conservation.) For reference, Fig. |l| shows the weak rates as a function of 
temperature. 
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2.1 Rate Expressions 

Generally, the weak-interaction rates may be expressed as integrals over phase space. For 
definiteness, consider the process ep — > vn. The other processes involve similar expressions. 
The rate, reactions per incident nucleon per time, is given by a twelve- dimensional integral, 

T ef ™ = - / II ^(2vr) 4 \Mf 5 (4) (e + P-V- n)f e f p (l - f v ){\ - f n ), (1) 
n p J i 

where n p is the number density of the incident nucleon, dUi = d 3 pi/[(27r) 3 2Ei] is the Lorentz 
invariant phase-space element, \A4\ 2 is the weak-interaction matrix element summed over 
initial and final spins (see Appendix A), e, p, v and n are the four- momenta of the particles, 
and the delta function 5^ (e + p — v — n) expresses conservation of four- momentum. In the 
rest frame of the thermal radiation, the phase-space densities f e , f p , f u and /„ are given by 
the usual Fermi-Dirac distribution functions. 

Assuming only the isotropy of space, we can integrate over dQ e , the direction of the 
electron's three-momentum and d(f) p , the azimuthal angle of the proton's three-momentum. 
After applying conservation of three-momentum to eliminate d 3 p n and conservation of energy 
to eliminate dp u the expression simplifies to the following five- dimensional integral: 



T ep ^ un = 29 ^ 6n J dp e dp p d cos 9 p d cos 9 u d<f) u 



x l¥lH \ M \ 2 f^ ~ /"X 1 " U ( 2 ) 

i±i e i±i p i±i n j 



(3) 



where E e , E p , E u , and E n denote the energies of the respective particles and J is the 
Jacobian introduced in integrating the energy part of the delta function. We take p e parallel 
to the polar axis (9 e = (f> e = 0) and (fi p = 0. We need an expression for E v in terms of the 
integration variables p e ,P P , @ P , @u, and <p u ; it is given by 



A 2 B + 2£ V M 4 - ml(AE 2 - B 2 ) 

Vv = 4E 2 - B 2 ' 

A 2 = 2E e E p + m 2 — m 2 —m 2 — m 2 — 2p e p p cos 9 P , 

B = 2 [p e cos 9 U + p p (cos 9 P cos 9 U + sin 9 P sin 9 U cos <fi u )] . (4) 

where E = E e + E p . If m v is taken to be zero (or very small) these equations simplify 
further. We left them in this form because for two-body processes in which the electron is 
the outgoing lepton, m u —^m e ^0. 
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2.2 Infinite-Mass Approximation 

If the nucleons are assumed to be infinitely massive,[] two things happen. The kinematics 
simplifies because the kinetic energy of the nucleons can be neglected (E e — E v = Q), and 
the matrix element simplifies, \M\ 2 = 2G 2 F (1 + 3g A )2E e 2E p 2E u 2E n . Here G F = 1.166 x 
10~ 5 GeV~ 2 is the Fermi constant and qa — 1.257 is the ratio of axial to vector coupling to 



the nucleon. The rate expression now reduces to the familiar form |31 

F oc G|(l + 3giK r 6(e 2 -g 2 ) 1 / 2 

ep " un 2tt 3 J q [l+exp(ez)] [1 + exp((g - e)z„))] ' [) 

where T is the photon temperature, T v is the neutrino temperature, e = E e /m e , z = m e /T, 
and z v = m e /T u . Expressions for the other five n p processes are similar. This is the 
approximation used in the standard nucleosynthesis code. 

In the era preceding nucleosynthesis, when to a good approximation baryons exist as free 
neutrons and protons, the neutron fraction X n is governed by 

X n = —X n T n ^ p + (1 — X n )T p ^ n , (6) 

where X p = 1 — X n , and 

r = r 4- r + r 

r?i— >p — r en^tvp ^ vn^tep ^n-^peu- (7) 

In thermal equilibrium, which holds when the rates are much greater than the expansion 
rate of the Universe, 

(-\ = Tp ^ n = ( ^\ e~ Q/T . (8) 
\pj T n ^ p \m p j 

In the infinite-mass limit, m n /m p = 1, and 

(n\°° _T^_ e _ Q/T ^ (g) 



■p) r 



oo 
n— »p 



The second equality in both equations is just detailed balance, which is one of the checks 
that we will use to test our results. 



3 Finite-Mass Corrections 

We carried out the five- dimensional integration using the Monte Carlo technique. The chal- 
lenge was to compute a finite- mass correction, whose size is of order a few percent, to a 
relative accuracy of a few percent. Hence, we had to evaluate the rate integral to order 

1 More precisely, m nj m p — > oo with m n — m p — Q fixed. 
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Figure 2: Finite-nucleon-mass correction to up — > en at T = 1.0 MeV as a function of 1/mjy. 
The error bars represent the estimated statistical uncertainty. The arrow indicates the actual 
value of the nucleon mass, ~ 940 MeV. Note, the correction vanishes for mjy — ► oo and 
is linear in 1/mjy. 

0.1% accuracy. The statistical uncertainty in integrating the rate expressions was typically 
around 0.3% for several minutes of computer time. Since this scales as 1/y/N, where N is 
the number of function evaluations, we could have achieved the accuracy requirements by 
simply increasing the integration time by a factor of ten or so. However, we found a more 
efficient way of proceeding, which also allowed other checks to be made. 

For E/rriN <C 1 the finite- mass correction to the weak rates should vary linearly with 
1/mjv (E ~ 1 MeV is the characteristic energy scale: E ~ Q,m e ,T.) As Ejm^ — > 0, the 
correction must approach zero. We adopted the following strategy. We performed a series 
of runs where we fixed the temperature but varied the nucleon mass. We found the finite- 
mass correction (= 5T = T — r°°) by applying a linear fit to the data and interpolating to 
rriN — 940 MeV. Figure shows the result of this procedure for the process ep — > un. The 
linearity of ST/T in 1/mjv is manifest as is the intercept at 6T /T = for — * oo. 

We found that the finite-mass correction for each process varies linearly with temperature. 
Therefore, we applied a linear fit to the correction as a function of temperature. In the end we 
were able to achieve the required order 0.1% absolute accuracy in the rate integrals. Figure |3] 
shows the finite-mass corrections <5r/r for all six processes that interconvert neutrons and 
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Figure 3: Finite-nucleon-mass corrections to the weak rates as a function of temperature, 
protons. The corrections are, as expected, of order a few percent. 



3.1 Tests of finite-mass corrections 

We applied several checks to our results. The simplest was to let 1/m/v vary. As 1/mjv — > 
the finite-mass rate corrections should approach zero; Fig. illustrates that they do. While 
Fig. illustrates convergence for a particular process and temperature, convergence was 
observed for all processes and temperatures. 

Another test is provided by detailed balance. Detailed balance requires that 

f^p— *n\ f n \ ( m n\ -Q IT ( ^ Q 



\r— = ; = Ur (10) 

V 1 n^p ) \P J e q \ m p / VZ nl N / 

while the neutrinos are coupled to the photon plasma and T v = T 7 . We numerically deter- 
mined r n _»p/r p _»„, and expressed it as 

r_^^ (3Q_ a + 1 \ e -Q/T. 



r n ^ p V 2 m N 

a = 1 corresponds to detailed balance being satisfied. Figure f| shows that detailed bal- 
ance is satisfied while the neutrinos are coupled to the photons and also when the neutrino 
temperature is set equal to the photon temperature. 
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Figure 4: Detailed-balance test. For a = 1 detailed balance is satisfied (see text). For the 
boxes, T v = T 7 is enforced; for the crosses neutrinos are assumed to be decoupled, so that 
T v < T 7 for T 7 < lMeV. Since the finite-mass corrections are obtained from fits to l/m^ 
and T 7 , the errors are strongly correlated. 

Finally, we were able to derive an independent and simpler expression for the finite-mass 
rates but using the infinite- mass matrix element. This provides a check of the complicated 
kinematics that arise for itln ^ oo. In Eq. (|l|) the three- momentum part of the delta function 
is expanded in complex exponentials, 

T ep ^ n = — (\{dU i £x{2n)e i ^ + ^-^-^ ) \M\ 2 

n p J i 

x5(E e +E P -E V - E n )f e f p (l - f v )(l - f n ). (12) 

If the integrand has no angular dependence, as is the case for the infinite-mass matrix 
element, the expression can be reduced to three-dimensional integral. (With angular de- 
pendence, the reduction can be done term by term in the angular expansion of the matrix 
element, but the resulting expression becomes very complicated.) For the infinite-mass ma- 
trix element all of the angular integrals can be done, 

, G 2 F (l + 3g 2 A ) frj , . . *d 3 x 
ej ™ = ^ST J [[dpiPiSm{x Pi ) — 

xd(E e + E P -E U - E n )f e f p (l - /„)(! - f n ). (13) 
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The integral over dx can be done by a standard contour integration, 

^lp-^un = F L 5 — I dE e dE p dE u E e E p E u (Ep + E e — E v 



bp 

-/-)(!-/»). ( 14 ) 



i=l 

where the energy delta function has been used to carry out the dE n integral and 

8 

= \Pe+Pp+Pu+Pn\ + \Pe+Pp-Pu-Pn\ 



i=l 



+ \Pe -Pp+Pu-Pn\ + \Pe ~ Pp ~ Pu + Pn\ 
-\Pe+Pp+Pu~Pn\- \Pe + Pp ~ Pu + Pn\ 

-\Pe-Pp+Pu+Pn\- \Pe ~ Pp ~ Pu ~ Pn\ ■ (15) 

Using standard techniques we carried out the numerical integration. 

For our test we compare T 3 with the rate obtained by inserting the infinite-mass matrix 
element into Eq. (0). This rate, denoted by T 5 , is given by 

F 5 ep -> vn = F 94 g — / dp e dp p d cos 6 p d cos 9 v d<j> v 

xplp 2 pPuE^f e fp(l - - f n ). (16) 

Formally, T 3 and T 5 are identical. Numerically, they are computed using independent tech- 
niques. Comparing them provides a stringent test of kinematics. Figure |5] shows this com- 
parison as a function of with T = 5 MeV. Figure |^ shows the comparison as a function 
of T with rriN = 100 MeV. In both cases, T 3 and T 5 agree within estimated numerical uncer- 
tainties. In conclusion, because our finite-mass rate corrections pass these three important 
tests, we have confidence that they are correct and accurate. 



4 Changes to the Standard Code 

The standard code calculates the weak rates at each temperature step. (The user is offered 
the choice of either integrating the rates of using fits to them. We chose to have the code 
integrate the rates at each step.) The rates are obtained by integrating equations like Eq. @, 
which assume that the nucleons are infinitely massive. Our approach was to implement the 
finite-mass rate corrections as a multiplicative factor at each temperature step: 




(17) 
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Figure 5: Kinematics test for T = 5.0MeV (see text). Circles are the five-dimensional 
integration and crosses are the three-dimensional integration. 



where 5T = J2i 8Tu r = J2i an d i runs over the three reactions that convert neutrons to 
protons (or vice versa). A plot of the relative finite-mass rate corrections as a function of 
temperature is shown in Fig. [7[ At the crucial freeze-out epoch (T ~ 0.8 MeV), both rates 
are corrected downward, with 5Y p ^ n /Y p ^ n = —0.43% and 5T n ^ p /T n ^ p = —0.55% 

One subtlety that we have not mentioned is the fact that in the standard code the 
weak rates are normalized to the measured free-neutron decay rate. This means that free- 
neutron decay has already been "corrected" for the zero-temperature, finite-nucleon-mass 
effect. Since the other rates are normalized to r n , they too include a piece of the finite- 
nucleon-mass correction, where the finite-nucleon-mass correction to free-neutron decay is 
5 n . To properly implement our results in the standard code we must "remove" this correction, 
so that Eq. (|17]) is actually implemented as 

r p _ = n", ( i + *pl (i - 5 n ) 



r 

5Y^ 



To compute S n we began with the twelve- dimensional phase-space integral for neutron 
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Figure 6: Kinematics test with = lOO.OMeV (see text). Circles are the five-dimensional 
integration and crosses are the three-dimensional integration. 



decay, 



IV^e, = - / WdYl^y \M\ 2 S^(e +p-v- n)f n (l - f p )(l - / e )(l - /„). (19) 



In the zero-temperature limit, f n = (n n /2)(27r) 3 5( 3 )(p n ), and (1— f p ) = (1— f e ) = (l—f v ) = 1. 
This expression can be simplified to the following two-dimensional integral: 



1 r Q-{Q 2 -ml)/2m N ,„/•!, „ \M\ p e E v 



n^peu 



E v 



2 7 7T 3 j me 

E v + p e cos 6 U 



i: 



dE P d cos 9, , 



m n E p \J\ 



1 + 



m n — E v — E e , 

m l ~ m l + m e _ 2m„E e 
2 (m n - E e +p e cos9 u ) 



(20) 

(21) 
(22) 
(23) 



where \J~\ is the Jacobian introduced due to integration over the energy delta function. Note 
that the full matrix element is used. Using standard techniques, we integrated this expression 
and found the zero-temperature neutron decay correction, 5 n = Y n ^ pev — = —0.206%. 
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Figure 7: Corrections to T p ^ n (T) and T n _, p (T). 

5 Results and Conclusions 

Our results, which were obtained for three massless-neutrino species and a mean neutron 
lifetime of 887 sec [32], are shown in Fig. |[ Plotted are the fractional changes in the mass 



abundances of D, 3 He, 4 He and 7 Li versus the baryon-to-photon ratio rj. The relative change 
in 4 He is approximately 0.47% - 0.50% over a large interval, 10~ n < rj < 10~ 8 . Over the 
most interesting part of the interval, rj = (2 - 7) x 10~ 10 , AY/Y = (0.49 - 0.50)%. Over 
the full interval in rj, the fractional change in the other light elements span the range 0.06% 
to 2.8%. Unlike 4 He, the inferred primordial abundances of these elements are known to 
nowhere near this accuracy and these changes are of little relevance at present. 

Several years ago the effect of finite nucleon mass on the yields of nucleosynthesis was esti- 



mated by Gyuk and Turner |33[ based on an approximation scheme developed by Seckel 
Seckel estimated the lowest-order correction, i.e., terms of order 1/mjv, to the weak rates 
and Gyuk and Turner incorporated them into the standard code. Our finite-nucleon-mass 
corrections are generally consistent with Seckel's estimates for the changes in the rates; how- 
ever, there are significant qualitative and quantitative differences. Likewise, our calculation 
of AY is consistent, but significantly smaller, than that of Gyuk and Turner who found 
AY ~ 0.006F based on Seckel's scheme. 
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Figure 8: Finite-mass corrections to light element abundances. For 4 He the fractional change 
is in the mass fraction; for the others it is in the number relative to H. 
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A Appendix: The Matrix Element 

Here we present the full matrix element to leading order in the weak interaction coupling 
constant |34|, |3"5|] . The particles labeled by 1 and 3 are the incident and outgoing leptons, 
and by 2 and 4 are the incident and outgoing nucleons. The matrix element is expressed in 
terms of the relativistic invariants s and t. Here, s — (1 + 2) 2 = (3 + 4) 2 and t — (1 — 3) 2 = 
(2 — 4) 2 , where 1, 2, 3, 4 are the four-momenta of the respective particles, ji = 1.81 is the 
anomalous weak magnetic moment of the nucleon, and g& = 1.257 is the ratio of axial to 
vector coupling of the nucleon. 

This expression applies to all six weak processes. For anti-leptonic, two-body processes, 
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e.g., e + n — * v e p, qa — > — Qa- For n — > pez/, s = (— 1 + 2) 2 and t = (— 1 — 3) 2 and for 
peis -> n, s = (1 + 2) 2 and t = (1 + 3) 2 . 

|A<| 2 = G|,(ti+t 2 + *3 + *4 + t5 + t6), 

*i = 16/ 2 /m^ x 

[—m\m 4 + m\m\ — m\m\m 4 + m 2 m 2 m 2 — m 2 m 2 s 

— m 2 m 2 t + m\m\m 4 + m\m 4 s — m\t + m\m 4 t 

— m 2 mg — m 2 m 2 m 2 + m 2 mgS + m 2 m 2 t + m 2 t 2 

+77137714 — m\m 4 s — m\m 4 t — m\t + m 4 t 2 ), 
< 2 = 8x 

(m\m\ — 2m 2 77i 2 77i4 + 2m\m\ + m\m 2 4 — 2m 2 s 
— m 2 t + 777 2 m 2 + 2m 2 m 2 — 27772S — m 2 t 
— 277i 2 m3m4 + 2m2m 4 t + m\m 2 4 — 2m\s — m 2 t — 2m 2 s 
-m 2 t + 2s 2 + 2st + t 2 ), 
t 3 = 4/|/m^ x 

(m\m\ — 2m\m,2m 4 — 3mfm 2 — mft + 2m 2 m 2 

+2m 2 m27n 2 — Am\m\s — 3m\m\t + Am\m2m' 3 m 4 — 2m 2 m 2 m 4 t 

+2777 2 77i 2 m 2 — 2m 2 mlt — 2m\m\ + Am\m 2 4 s + m\m 2 4 t 

+Am\st + m 2 t 2 — 2m2777 2 — 2m 2 £ — Sm^wij 

+4m 2 777 2 s + m\m\t + mm\st + 2m2t 2 — 2m2m 3 m 4 

0/10/1 o 

—2rri2m 3 m 4 t + Am2m 4 t + m 3 m 4 — m 3 t + Am 3 st 
+m 2 3 t 2 - 2m\t + 4m 2 st + 2m 2 t 2 - 4s 2 t 
-4 S t 2 ), 

f 4 = %1 X 

(m 2 m 2 + 2m 2 m 2 m 4 + 2m\m\ + m\m 2 4 — 2m\s 
—m\t + 7n 2 m 2 + 2m 2 777 2 — 27772S — 777 2 £ 
+2m27n 2m 4 ~~ 2m.2m 4 t + m 2 m 2 — 2m 2 s — 7n 2 t 
-2m 2 s - m 2 t + 2s 2 + 2st + t 2 ), 

/ 23 22 1 2 2 1 2 j. , 2 3 

( — m^m^ — m l m 2 m 4 + m 1 m 2 m 4 + m 1 m2i + m x m\ 
+m\m 4 t + m\m\ + m^t + 777.2777.3777.4 + mlm 4 t 

— 7772777 2 777 2 + 7772777 2 £ + 7772777 2 t — 2777 2 St — 777 2 t 2 
— 777 3 777 4 + m 3 m 4 t + 777 4 t — 2m 4 St — m 4 t ), 
t 6 = X 
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{—m\m\ + m\m\ + m\t + m\m\ + m\t 
—m\m\ + m\t + m^t — 2st — t 2 ). 
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